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tion the evolution of a system is, in general, extremely sensitive to initial conditions 
and to small changes in system parameters. This sensitivity is an intrinsic physi- 
cal property of the system, and not merely of numerical solutions. Sensitivity to 
initial conditions - aka chaos - may be found in the conservative system also, but 
not necessarily in all regions of phase space. In contrast, in the presence of weak 
dissipation, a nearly-Hamiltonian system is forced to sweep across large regions of 
phase space, including chaotic regions near resonances. Consequently, sensitivity of 
the evolution is always present in this case, for any initial conditions. The scatter 
in Figure 3 is indirect evidence of this. It is related to the complicated dynamics 
of resonance passage phenomena. A discussion of this phenomenon is beyond the 
scope of this paper. However, it is worth pointing out that resonance passage can 
produce "anomalous" behavior, i.e. evolution that is characterized by singularly 
large changes of orbital parameters in relatively short times. This phenomenon 
is well-known in the literature and is an active area of current research (see, for 
example, Malhotra 1994, for a recent review). The meaning of numerical solutions 
in the presence of such dynamical behavior is also an interesting question (see, for 
example, Quinlan & Tremaine 1992). 
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4. SUMMARY 

In this paper, we have introduced a new numerical integration method for Solar 
system problems with small dissipation. This method is a modification of the 
second-order "N-body map" introduced recently by Wisdom & Holman (1991) 
(also independently by Kinoshita et al. (1991)) for nearly-integrable Hamiltonian 
systems of the form H = -ffkepler + I^H' where -ffkepler and /^H' are separately 
integrable. -ffkepler describes the unperturbed Keplerian motion and /^H' describes 
the small mutual gravitational perturbations of the planets (or satellites, in a 
planetary satellite system). 

The new method is summarized as follows. Recall that for the conservative 
problem, a single step of the N-body map with size r requires three operations. 
These are: (i) evolve the system according to -ffkepler for time ^r; (ii) evolve the 
system according to /^H' for time r; (iii) evolve the system according to -ffkepler 
again for time ^r. The evolution of the unperturbed system, according to ifkeplen 
is efficiently accomplished with the use of Gauss' / and g functions to update 
the position and velocity of a particle on a Keplerian ellipse. Because the mutual 
perturbations described by /^H' are a function of positions alone, sub-step (ii) 
involves only small increments ("kicks") in the velocities. In the presence of a 
small dissipative force, Ar = eF, we have shown how this integration scheme can 
be used with modified functions / and g, and a new function, h, to update the 
positions and velocities in sub-steps (i) and (iii). The corrections, f — f, g — g and 
h are given in terms of the drag force, eF. 

The method we have described here is a second order integrator: the leading 
truncation error in the positions per time step r is 0{iit^) + 0{£t^). We have 
verified this numerically by applying this new method to the restricted three-body 
problem with a small gas drag force on the test particle. We have given a com- 
parison of the numerical errors of this method with the fourth order Runge-Kutta 
integrator. The numerical tests show that the position errors grow linearly with the 
length of integration (in contrast with the quadratic rate of error growth with the 
Runge-Kutta integrator). Furthermore, the spurious "numerical dissipation" that 
is present in conventional integrators such as the Runge-Kutta is nearly absent in 
the 'mapping', as shown by the error in the Jacobi parameter. From the results of 
numerical tests, we estimate that the 'mapping' is at least an order of magnitude 
faster than the fourth-order Runge-Kutta for similar levels of integration error. 

An obvious disadvantage of the "generalized mixed variable mapping" described 
in this paper is its low order. For conservative systems, symplectic integrators 
of higher order are described in Yoshida (1990), Wisdom & Holman (1991) and 
Kinoshita et al. (1991), and other refinements that yield significantly improved 
performance are described in Saha & Tremaine (1992). The feasibility of similar 
improvements in the generalized mapping for dissipative systems will be considered 
in the future. 

Finally, we caution the reader that in planetary dynamics problems with dissipa- 
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Fig. 3. The inaxiinuin error in the Jacobi integral, C, and in the mean longitude, A, after an 
integration for 10,000P are shown as a function of the gas drag parameter, K. These integrations 
were done with the 'mapping' with step size r = P/100. For reference, a linear relationship 
between the abcissa and the ordinate is indicated by the dashed lines. 



that |5A|inax and |5C|niax both increase as ~ r"*'^ for RK4.^ However, the numerical 
factor is dramatically larger than for the 'mapping'. Consequently, the 'mapping' 
method easily out-performs RK4 for larger step sizes. 

In Figure 3 we show the maximum integration error as a function of the magni- 
tude of the drag coefficient, K . There is considerable scatter in this plot, but the 
trend is a roughly linear increase of the error with K. This is as expected because 
in our integrations, the (9(er"^) positional error from the drag force is generally 
comparable to or exceeds the 0{iit^) error from the planetary perturbation. (In 
our numerical tests ji ~ 10"'*; the value of e (cf. Eqn. 18) varies considerably 
during each run, but its maximum value is approximately 0.01/i'.) 



It is surprising to find an exponent 4.5, rather than 4 for a fourth-order integrator. We do 
not have an explanation for this better-than-expected performance of RK4, but see the discussion 
in Press et al. 1989. 
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Fig. 2. The inaxiinuin error in the Jacobi integral, C, and in the mean longitude, A, after 10,000 
planet orbits when the orbit of Fig. 1 is integrated with different step sizes, r. (a) and (b) show 
the errors incurred with the 'mapping'; (c) and (d) show the errors incurred with a fourth order 
Runge-Kutta integrator. 



the 'mapping', and the panels (c) and (d) on the right are for the Runge-Kutta 
integrator. In these integrations we set the drag coefficient K = 0.01. 

In Figure f we plot the error in C and in the mean longitude, A, against time 
for a run using step size r = P/fOO. We observe that the mean longitude error 
grows linearly with time in the 'mapping' method, but quadratically with time in 
the RK4 integration. Furthermore, with the 'mapping', there is virtually no secular 
growth of the error in C (recall that this is also characteristic of the mixed- variable 
symplectic integrators for the conservative system); with the RK4 integrator, the 
error in C grows linearly with time. This shows that spurious "numerical dissipa- 
tion" is virtually absent in the 'mapping'. 

Figure 2 shows the maximum errors incurred as a function of step size. We 
observe that with the 'mapping', the position error increases as ~ r^ as does the 
error in C. This shows that the 'mapping' method is of second order. Note also 



computer time as one step of the RK4. Thus, for the same step size, the 'mapping' is approximately 
twice as fast as RK4. 
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Fig. 1. A test particle orbit subject to a planet's gravitational perturbations as well as gas drag 
forces is integrated for 10,000P, where P is the orbital period of the planet. (The initial conditions 
and other parameters adopted in this integration are given in the text.) The absolute value of 
the errors in the instantaneous Jacobi integral, C, and in the instantaneous mean longitude, A, 
are plotted against time (indicated in units of P). (a) and (b) show the errors incurred with the 
generalized 'mixed variable mapping' described in the text; (c) and (d) show the errors incurred 
with a fourth order Runge-Kutta integrator. These integrations were performed with a fixed step 
size T = P/100. 



done for different values of K from 0.0 to 0.1. To determine the errors in an 
individual orbit integration, we ran a second integration using a smaller step size. 
The smallest step-size, r, used was 1/1600 of the planet's orbital period, P. The 
errors displayed in Figures 1-3 compare the numerical solution using the indicated 
step size with a numerical solution obtained with r = P/1600. 

In Figures 1 and 2 we display a comparison of the 'mapping' with a fourth 
order Runge-Kutta (RK4) integrator.^ The panels (a) and (5) on the left are for 



^ We note that RK4 makes four force evaluations per time step, whereas the 'mapping' makes, 
formally, only one force evaluation per step (and this does not include the central force). However, 
moving the particle with the /, g and h functions consumes slightly more machine time than one 
force evaluation in RK4. (This estimate can be expected to be somewhat sensitive to the form of 
the dissipative force.) We find that one time step with the 'mapping' uses approximately half the 
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where Mq is the mass of the Sun, and r, v are the heliocentric position and velocity 
of the test particle; the heliocentric position of the planet, Fp is given by 

rp(i) = ap(cosrapi, sinrapi), (15) 

and nip, Up and Wp are the planet's mass, orbital radius and orbital frequency, 
respectively. The gas drag force per unit mass acting on the particle is taken to 
be that appropriate for high Reynolds number flow (cf. Landau & Lifshitz 1976, 
Adachi et al. 1976): 

eF = -Kuu, u = V - Vgas, (16) 

where Vgas is the velocity of the gas which is assumed to be in Keplerian rotation, 
but partially supported by internal pressure: 

Vgas^ (l-??)vkep(r). (17) 

Here Vkep(r) is the local Keplerian velocity, Ukep = ^GM^g/r, and 77 is a measure 
of the radial pressure gradient in the gas. In models of the Solar Nebula, r] is 
a few times 10""^ (cf. Weidenschilling & Cuzzi 1993). The drag coefficient K is 
proportional to the density of the gas and to the surface-area-to-mass ratio of the 
particle; it has dimensions of inverse length. It is worth pointing out that the ratio 
of the drag force to the central force is given by 



e = rK[ . (18) 

For a circular orbit, w/ukep ~ 1]; for orbital eccentricity e exceeding r], w/ukep ~ e. 
Therefore, in most cases of interest, (w/ukep)^ <€i 1, so that even fairly large values 
of K represent a small perturbation to the underlying Keplerian orbit. 
For K = 0, there exists a constant of the motion, the Jacobi integral, 

C=-/7 + np(rx v),. (19) 

The value of C determines the global stability properties of the test particle orbit. 
It is therefore a natural parameter - analogous to the total energy and angular 
momentum for the general N-body problem - for monitoring the error of a numer- 
ical solution. For K / 0, C is not constant; nevertheless, we can define a value of 
this parameter at any instant of time during the evolution of test particle. In the 
numerical test described below, we quote the errors in C and in the mean longitude 
of the test particle. 

We choose units such that G = Mq -|- nip = Op = 1. In these units, the planet's 
orbital frequency is unity and its orbital period is 27r. In the numerical tests, the 
initial conditions of the test particle were a = 1.5, e = 0.1, uj = 7r/2, X = {uj 
and A being the longitude of perihelion and the mean longitude, respectively). We 
used nip = 10"'*, and started the planet at conjunction with the test particle at 
zero longitude. The gas drag parameter r] was fixed at 0.005. Integrations were 
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Substituting Eqn. 10 in Eqn. 9, we obtain 

/ 1 oGM 1 .aCMro-vo \ / 1 .GM x 



+e 



r Q U / Q ' U / Q 



2 ^^10 + g^^io + • • • jro + [-t'F2o + -r-'i^ao + • • • )vo 



+ (2 ^^^30 + qt'^Fso + • • • ) (ro X vo) 



(12) 



Comparing with Eqn. 6, it is easy to see that, as expected, the (9(e°) coefficients of 
ro and vq are the Taylor series expansions of Gauss' / and g functions. And, finally, 
to leading order in r we have the following expressions for the 0(e) corrections: 

f{to + r, to) - f{to + r, to) = -t\Fw, 

g{to + T,to) - g{to + T,to) = -r^eF2o, (13) 

h{to + T,to) = -T^eFso- 

Provided F{t,r,v) is a sufficiently smooth function, the terms higher order in r 
can be made negligibly small with suitably small step sizes. 

Straightforward bookkeeping shows that the leading order truncation errors 
per time step with the above scheme are 0{iit^) + 0{£t^) in the position, and 
0{iit^) + 0{£t'^) in the velocity, where the first term derives from the gravitational 
perturbation jiH' and the second from the drag force. Since the position error per 
time step is (9(r"^), the order of the integration is formally 2. 

3. NUMERICAL RESULTS 

In this section, we apply the new method to one of the simplest model problems 
in Solar System dynamics with dissipation. We give a detailed numerical error 
analysis, comparing the new method with the standard fourth order Runge-Kutta 
method (RK4). 

The model problem we have chosen is the planar circular restricted three body 
problem with gas drag. The system consists of a test particle in a heliocentric orbit 
perturbed by a planet; the test particle orbit decays under the infiuence of a gas 
drag force whose magnitude is proportional to the square of the relative velocity 
of the particle with respect to the ambient gas. The gravitational perturbations 
on the particle are described by the Hamiltonian 

1 o GMr. ^ / 1 r-rp\ 

2 r Vr — Fr, rt 
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Towards this end, we note that for an arbitrary eF, the motion, although nearly 
Keplerian, need not remain in a plane. Therefore, we introduce a generalization of 
Eqns. 2 as follows: 

(6) 
v(0 = f{t,to)ro + g{t,to)vo + h{t,to){roXvo), 

where the modified f and g functions and the new function, h, take account of the 
external dissipation. One can see directly that f — f, g — g and h are all 0(e). 
Below we derive expressions for these 0(e) corrections. 

Without loss of generality, the vector F can be written as follows: 

F(i, r, v) = Fi (i, r, v)ro + F2 (t, r, v) vq + F3(t, r, v) (fq X vq) , (7) 

where 

1 rro -F - (ro • vqIvq -Ft 



Fr(M,v) = -[ 



F3(t,r,v) 



ro L 1 - (ro ■ vq)"^ 

1_ r vo -F - (ro --^0) ^0 • F 
uo L 1 - (ro • vo) 

(ro X vo) • F 
|ro X vop 



where ro and vo are the unit vectors along ro and vo, respectively. 
Consider the Taylor series, 

r(to + r) = ro + rro + -r^ro + -t^ r + • • • (9) 

Z D 



Now, 



ro = Vo, 

GM 

'■o = — ;^ 



r, 



ro + e[Fioro + F20V0 + F3o(ro X vo)] , (10) 





Vo 3(ro-vo)ro 



^n^r'^o ^iro-vojro] , f/, , /, , /^ / ^ 

ro = -GM\— ^ 5-^ +eFioro + F2oVo + F3o(ro X vo, 

LTq ''o 

where 

Fko = Fk(to,ro,vo), 



d d d 

Fko = ^Ffc(io,ro,vo) + ro • ^-Ffc(io,ro,vo) +ro • ^--Pfc(io,ro,vo) 
ot or av 
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can intuitively recognize the advantage of this method: the unperturbed Keplerian 
motion is not completely reinvented at each step, as is the case in conventional 
integrators. 

The algorithm described above is for a second order integrator; the position error 
per time step is 0{iit^). (We follow the convention that an order n integrator has 
truncation error (9(r""'"^) per time step). The method can be generalized to higher 
order, and certain other refinements can be made. These improvements as well as 
the numerical stability and truncation errors of the MVS integrators are analyzed 
in detail by Yoshida (1990), Wisdom & Holman (1992) and Saha & Tremaine 
(1992). 

Having described the 'mixed variable symplectic' method, we now turn to a 
straightforward exploitation of this scheme for weak dissipative forces superposed 
on nearly-Keplerian motion. Consider the following equations of motion: 

. _ dH_ 

(3) 
dH 

p = -ir + ^^' 

or 

where H is as given in Eqn. 1, and 

eF = eF(i,r,v), (4) 

describes a small, "smooth", local dissipative force. 

The modification of the MVS method for the above system is as follows. In 
order to solve the initial value problem with (r, v) = (fq, vq) at time t^ to obtain 
(r(i), v(i)) at time i = ig + t, where r is small relative to the unperturbed orbital 
period, we retain the separation of H into -ffkepler and j-iH' ^ and also retain the 
"drift" -"kick" -"drift"... sequence in the integrator. The "kicks" to the momenta 
are the same as before in the conservative case (arising from the mutual gravita- 
tional interactions), but the "drift" part is no longer a purely Keplerian motion. 
That is, in-between the "kicks" we must now advance the positions and momenta 
according to -ffkepler P^us the effects of the perturbation, eF. Thus, we need to solve 
the following equations of motion in the first and third sub-steps: 



(5) 

GM 
V = —Y + eY 

where —GM/r represents the unperturbed Newtonian potential for the j'-th body. 
Accordingly, we seek a generalization of Eqn. 2 that accomplishes the "drift" be- 
tween the "kicks" to some desired accuracy. 
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studies involving long integration times. Therefore, computationally efficient meth- 
ods of integration for a few-body 'planetary system' with dissipation are always 
desirable. This reason has motivated the present work. 

In Section 2, we describe a straightforward modification of the N-body map that 
allows the inclusion of small non-Hamiltonian perturbations. The new method may 
be called a "generalized mixed variable mapping". (In the rest of this paper, we 
refer to this method simply as the "mapping" when there is no risk of ambiguity.) 
In Section 3 we describe an application of this method to the restricted 3-body 
problem with gas drag. We give detailed numerical error analysis of this problem, 
comparing the new method with the commonly used fourth order Runge-Kutta 
integrator for first-order ODE's (see, for example. Press et al. 1989). This study 
demonstrates the superior performance of the new method for planetary dynamics 
problems with weak dissipation. Section 4 is a summary. 

2. GENERALIZED MIXED VARIABLE MAPPING 

First, we give a brief description of the 'mixed variable symplectic' method in 
the implementation introduced by Wisdom & Holman (1991). Operationally, this 
technique may be summarized as follows. In order to integrate the system of Eqn. 1 
from time to to time to-\-T (where r is a small increment, typically less than 1/10 
of the smallest orbital period): 

— advance r, p according to -ffkepler for time ^r; 

— apply a "kick" of magnitude Ap = —t/^ (dH'/dr) to the momenta, p, and 

— advance the new positions and momenta according to -ffkepler for time ^r once 
again. 

We have dropped the index j for clarity. Note that the first and last sub-steps 
involve evolution on the unperturbed orbit (a Keplerian "drift"). This can be 
trivially accomplished in terms of the Keplerian orbital elements (or more generally, 
action-angle variables), while the second sub-step is trivial in cartesian variables. 
Hence the terminology, 'mixed variable'. As pointed out by Wisdom & Holman 
(1991), a convenient and efficient scheme to advance the positions and velocities, r 
and V = p/m, in-between the "kicks" (when the motion is strictly on a Keplerian 
ellipse) is encapsulated in Gauss' / and g functions (see Battin 1987, or Danby 
1988): 

r(0 = f{t,to)ro + g{t,to)vo, 

(2) 
v(0 = f{t,to)ro + g{t,to)vo. 

This solution is possible because (r, v) and (fq, vq) are non-collinear vectors that 
lie in the same plane. The / and g functions require a determination of only a 
few of the osculating orbital elements from the initial conditions (ro,vo). One 



2 RENU MALHOTRA 

and position errors grow quadratically with time. For this simple reason, symplec- 
tic integrators have a significant advantage over conventional integrators for long 
integrations of a planetary system. 

Since the first application by Wisdom & Holman (1991) to integrate the outer 
planets for a billion years (obtaining results that compared very well with the 
earlier results of the 'Digital Orrery' calculations of Sussman & Wisdom (1988)), 
this method has been used for several other problems in Solar System dynamics. 
Sussman & Wisdom (1992) followed all the planets for 100 million years; Saha 
(1992) adopted this method for a numerical investigation of the chaotic dynamics 
of asteroids near the 3:1 Kirkwood Gap; Levison & Duncan (1993) and Holman 
& Wisdom (1993) have used it to explore the long term stability of test particle 
orbits in the outer Solar System. 

While the method is demonstrably very useful for a certain class of problems 
in Solar System dynamics, it is obvious that the Hamiltonian system of Eqn. 1 
is a highly idealized description of an N-body system — that of point masses 
subject only to mutual gravitational interactions according to Newtonian gravity. 
Although the planetary system today resembles closely this idealized Hamiltonian 
system, this was not the case in the early stages of the Solar System. Even if 
one confines oneself to the present state, for the dynamical behavior over billion 
year timescales one may worry about additional physics (cf. Nobili et al. 1989) — 
the effects of higher order gravitational moments of the bodies, general relativistic 
effects, mass loss of the Sun by radiation and the solar wind, etc. In questions 
related to the early stages of the Solar System and planet formation, point mass 
interactions provide the dominant physics, but other dissipative processes such 
as inelastic collisions and gas drag are also important. Moreover, there is a large 
class of problems concerning the contemporary Solar System, such as the ring 
systems and satellite systems of the planets, in which forces other than point mass 
gravity are quite important even when relatively small. Here the quadrupole and 
higher order moments of the central planet, and collisional or tidal dissipation have 
important qualitative effects on the long term evolution. Yet another example is 
cometary orbital evolution where non-gravitational effects such as outgassing are 
small but of great interest. As a final example, we can mention the evolution of very 
small particles (dust grains) in the Solar System which are infiuenced significantly 
by radiation forces. 

Of the examples cited above, higher order gravitational moments are straight- 
forwardly incorporated in the 'mixed variable symplectic' schemes because (i) they 
do not change the symplectic nature of the system, and (ii) they still allow a sep- 
aration of the Hamiltonian into an unperturbed (Keplerian) and a perturbation 
part, each of which is separately integrable. However, dissipative forces cannot 
be accommodated within the 'mixed variable symplectic' schemes, for, in general, 
they destroy the symplectic structure of the equations of motion. Because weak 
dissipation is of interest in many Solar System problems, and many of those prob- 
lems are not amenable to analytical methods of study, this necessitates numerical 
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1. INTRODUCTION 



*~~5 Efforts to solve the question of the long-term stability of planetary orbits have 

l/~J recently spun-off a new algorithm (Wisdom & Holman 1991, Kinoshita et al. 1991) 

^""^ for the numerical integration of the conservative gravitational few-body problem 

described by a Hamiltonian of the following type: 

g H{rj, pj) = /7kepler(rj, Pj) + fiH'{rj). (1) 

r Here ifkepier(i'ij Pj) is the usual Hamiltonian for the two-body problem which de- 

~^ scribes the unperturbed Keplerian orbit of the j'-th planet, and fj,H'(rj) describes 

Q^^ the much smaller mutual gravitational interactions of the planets. The new tech- 

"^ nique exploits two features of the system: (i) the symplectic structure of the equa- 

O tions of motion, and (ii) the integrability of both the unperturbed Hamiltonian, 

' ifkepien and the perturbation Hamiltonian, i^H' , when taken separately. The latter 

5^ property permits an efficient algorithm based upon the use of two sets of canonical 

c/2 variables (one for -ffkepler and the other for jiH'). Kinoshita et al. call this method 

^ a "modified symplectic integrator", while Wisdom & Holman call it an "N-body 

map". Saha & Tremaine (1992) have reviewed a class of symplectic integrators 

that are particularly suitable for Solar System dynamics; they refer to these by 

the more descriptive name, "mixed variable symplectic" [MVS] integrators. 

The defining feature of the MVS integrators is that they preserve the symplec- 
ticity of the system and they do not produce a secular error in its first integrals. 
Therefore, quantities such as total energy and angular momentum (or the Ja- 
cobi integral for the restricted three-body problem) do not exhibit secular changes 
and the position error grows only linearly with time in a numerical solution ob- 
tained with these integrators. In contrast, conventional integrators such as the 
Runge-Kutta introduce a so-called "numerical dissipation" such that the numeri- 
cal solution exhibits a secular change of total energy even in a conservative system. 



